
%% Finite Slit Diffraction

lambda=4; %Wavelength
k=2*pi/lambda
a=20; %Width of Slit

%Simulation Area Size
Xmax=200
Ymax=50

%Use 'N' sources to represent the slit
N=100;

Svec=linspace(-a/2,a/2,N)
WaveAmp=0;

%Sum over point sources
for ss=1:N
    [X,Y]=meshgrid((0:.1:Xmax),(-Ymax+Svec(ss):.1:Ymax+Svec(ss)));
    R=sqrt(X.^2+Y.^2);
    Wave=exp(i*(k*R))./R; %From Huygens Principle
    WaveAmp=WaveAmp+Wave;
end


ImageScale = [-.1,.1];
%ImageScale = [];

figure
imshow(real(WaveAmp),ImageScale)

%% Time Average Intensity

lambda=4; %Wavelength
k=2*pi/lambda
a=20; %Width of Slit

%Use 'N' sources to represent the slit
N=10;
Svec=linspace(-a/2,a/2,N)

%For Time Average
w=1;
Tvec=linspace(0,2*pi/w,20);

WaveInt=zeros(size(R));

for tt=1:length(Tvec)
    t=Tvec(tt);
    WaveAmp=0;
    for ss=1:N
        [X,Y]=meshgrid((0:.1:Xmax),(-Ymax+Svec(ss):.1:Ymax+Svec(ss)));
        R=sqrt(X.^2+Y.^2);
        Wave=exp(i*(k*R-w*t))./R; %From Huygens Principle
        Wave=exp(i*(k*R-w*t)); %Approximation for Large distances
        WaveAmp=WaveAmp+Wave;
    end
    WaveInt=WaveInt+abs(WaveAmp).^2;
end

WaveIntAve=WaveInt/length(Tvec);

ImageScale = [];
%ImageScale = [-1,1];

figure
imshow((WaveIntAve),ImageScale)


%% Single Slit Far from Aperature

% Laser light, which can be modeled as a plane wave, is incident on a
% finite width slit.  The wavelength of light is 633 nm.  The slit
% width is 20 microns wide and the center of the slits are 100 microns apart.
% We use Huygen's principle to find the intensity of the light
% on a screen 10 cm away. 

%Use microns as units
Lambda=.633;
a=20;
L=1e5;

%Make our own sinc function
mysinc=@(x)sin(pi*x)./(pi*x);


%Setup an array of sources at opening of slits
Nsources=50;
YSourceVec=linspace(-a/2,a/2,Nsources)';

%Screen coordinates
YScreen=(-1e4:.1:1e4)';

%Function to calculate amplitude
wave_amp=@(YSource,Phase) 1./sqrt(L.^2+(YScreen-YSource).^2) .* ...
    exp(i*(2*pi/Lambda*sqrt(L.^2+(YScreen-YSource).^2)+Phase));

WaveInt=0
for Phase=(.1:.1:2*pi); %Time Average
    WaveAmpAll=0;
    for ii=1:length(YSourceVec) %Sum over sources
        WaveAmpAll=WaveAmpAll+wave_amp(YSourceVec(ii),Phase);
    end
    WaveInt=WaveInt+abs(WaveAmpAll).^2;
end
%normalize
WaveInt=WaveInt/max(WaveInt);

%Plot
figure
plot(YScreen,WaveInt)

%Theory
SinTheta=YScreen./sqrt(L.^2+YScreen.^2);
SingleSlitDif=mysinc(SinTheta/Lambda*a).^2;
hold on
plot(YScreen,SingleSlitDif,'--','linewidth',2)

xlabel('Screen Position (cm)')
ylabel('Normlized Intensity')
legend('Huygens'' Principle','Fraunhofer Diffraction Theory')


%% Double Slit Far from Aperture

%Make our own sinc function
mysinc=@(x)sin(pi*x)./(pi*x);

%Use microns as units
Lambda=.633;
a=20;
d=100;
L=1e5;

%Setup an array of sources at opening of slits
Nsources=50;
YSourceVec1=linspace(d/2-a/2,d/2+a/2,Nsources)';
YSourceVec2=-linspace(d/2-a/2,d/2+a/2,Nsources)';
YSourceVec=cat(1,YSourceVec1,YSourceVec2);

%Screen coordinates
YScreen=(-1e4:.1:1e4)';

%Function to calculate amplitude

wave_amp=@(YSource,Phase) 1./sqrt(L.^2+(YScreen-YSource).^2) .* ...
    exp(i*(2*pi/Lambda*sqrt(L.^2+(YScreen-YSource).^2)+Phase));


WaveInt=0
for Phase=(.1:.1:2*pi); %Time Average
    WaveAmpAll=0;
    for ii=1:length(YSourceVec) %Sum over sources
        WaveAmpAll=WaveAmpAll+wave_amp(YSourceVec(ii),Phase);
    end
    WaveInt=WaveInt+abs(WaveAmpAll).^2;
end
%normalize
WaveInt=WaveInt/max(WaveInt);

%Plot
figure
plot(YScreen,WaveInt)

%Theory
SinTheta=YScreen./sqrt(L.^2+YScreen.^2);
doubleSlitDif=cos(pi/Lambda*SinTheta*d).^2.*mysinc(SinTheta/Lambda*a).^2;
hold on
plot(YScreen,doubleSlitDif,'--','linewidth',2)

xlabel('Screen Position (cm)')
ylabel('Normlized Intensity')
legend('Huygens'' Principle','Fraunhofer Diffraction Theory')
